Internal models for nonlinear output agreement and optimal flow 

control 



Mathias Burger* 1 and Claudio De Persist 2 

institute for Systems Theory and Automatic Control, University of Stuttgart, 
| Pfaffenwaldring 9, 70550 Stuttgart, Germany 

q ■ mathias.buerger@ist.uni-stuttgart.de 
{Nj . 2 ITM, Faculty of Mathematics and Natural Sciences, University of Groningen, 

! Nijenborgh 4, 9747 AG Groningen, the Netherlands c.de.persis@rug.nl 

February 5, 2013 



- 1—1 

X 



Abstract 



CO 

^2 ■ This paper studies the problem of output agreement in networks of nonlinear dynamical systems 

under time- varying disturbances. Necessary and sufficient conditions for output agreement are de- 
rived for the class of incrementally passive systems. Following this, it is shown that the optimal 
distribution problem in dynamic inventory systems with time-varying supply and demand can be 
cast as a special version of the output agreement problem. We show in particular that the time- 
varying optimal distribution problem can be solved by applying an internal model controller to the 
00 ' dual variables of a certain convex network optimization problem. 

t> 
o 

csj ■ 1 Introduction 

O " 

' Output agreement has evolved as one of the most important control objectives in cooperative control. 

It has been studied in various different contexts, ranging from distributed optimization ([TBA86]) up 
to oscillator synchronization ([SS07]). Adding up to these results, we discuss in this paper the output 
agreement problem in the context of optimal distribution control for inventory networks with time- 
varying supply. 

Internal model control tools have been used to handle output agreement problems in a variety of different 
formulations, see e.g.,([WSAll]), ([BAW11]), ([DJ12b]). We consider here a different problem set-up, 
involving time- varying external disturbance signals, and solve the output agreement problem for the class 
of incrementally passive systems. Our derivations follow the trail opened by [PM08]. 
The output agreement problem with time- varying external signals, studied in this paper, turns out to be of 
particular relevance for the routing control in inventory systems. We consider a simple inventory system, 
taking into account the storage levels and routing between the different inventories. This dynamics 
models, e.g., supply chains ([AGT11]) or data networks ([MS83]). A key challenge in such inventory 
systems is to handle a time- varying supply and demand in an optimal way, while using only a distributed 
control strategy. 

The contributions of this paper are twofold. First, we present necessary and sufficient conditions for the 
output agreement problem under time-varying disturbances. We consider networks of nonlinear systems 
interacting according to an undirected network topology. Following an internal model control approach, 
we consider controllers placed on the edges of the network and provide necessary conditions for the output 
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agreement problem to be feasible. Sufficient conditions for output agreement in networks of incrementally 
passive dynamical systems are provided. We prove that the output agreement problem is feasible if one 
can find an incrementally passive internal model controller. As a second contribution, we show that 
the optimal distribution problem in inventory systems with time-varying supply and demand can be cast 
as an output agreement problem. The necessary conditions for the optimal distribution problem are 
a specific representation of the well known regulator equations. Subsequently, we present controllers 
solving the optimal distribution problem, for either quadratic cost functions or constant supplies. The 
internal model approach generalizes the results of [De 13]. 

The remainder of the paper is organized as follows. The basic problem formulation and necessary 
conditions for output agreement are presented in Section 2. Sufficient conditions for output agreement 
in networks of incrementally passive systems are discussed in Section 3. The time- varying optimal 
distribution problem is introduced in Section 4, where also necessary conditions are discussed. We 
present then the solution to the problem for linear-quadratic problems in Section 4.1 and for constant 
supplies in Section 4.2. 

Preliminaries: The notation we employ is standard. The set of (positive) real numbers is denoted 
by K (M>). The distance of a point q from a set A is defined as dist^g = inf pG _4 ||p — q\\. The range-space 
and null-space of a matrix B are denoted by 1Z(B) and M(B), respectively. A graph Q = (V,E) is an 
object consisting of a finite set of nodes, \V\ = n, and edges, \E\ = m. The incidence matrix B e M™ xm 
of the graph Q with arbitrary orientation, is a {0, ±1} matrix with [B]ik having value '+1' if node i is 
the initial node of edge k, '-1' if it is the terminal node, and '0' otherwise. We refer sometimes to the 
flow space of Q as the null space J\f{B T ), and the cut space of Q as the range space 1Z(B). Additionally, 
M{B) is named the circulation space of Q, and 1Z(B T ) the differential space. 



2 Problem formulation and necessary conditions 

We consider a network of dynamical systems defined on a connected, undirected graph Q = (V, E). Each 
node represents a nonlinear system 

Xi — fiiXii^ii^i) 

Vi = hi(xi,Wi), i = 1,2, ...,n, 

where Xi £ R n is the state, Ui,yi € W are the input and the output, respectively. Each system (1) 
is driven by the time- varying signal Wi <E representing, e.g., a disturbance. We assume that the 
exogenous signal Wi is generated by a dynamical system of the form Wi = Si(uii). In what follows it will 
be useful to require a gradient structure, i.e., Si(wi) = VEi(wi), with Y>i{wi) a concave function. Then 
for all Wi , w[ 

(wi - w-) T (s l (w 4 ) - Si(w'i)) < 0. 
As an example of VEi(wi) consider the linear function with skew-symmetric matrix 

VEi(tUi) = S iWi , Sf + Si = 0. 

We stack together the signals Wi, for i = 1, 2, . . . , n, and obtain the vector w, which satisfies the equation 
w = s(w). In what follows, whenever we refer to the solutions of w = s(w), we assume that the initial 
condition is chosen in a compact set W = W\ x . . . x W n . Similarly, let x, u, and y, be the stacked 
vectors of Xi 7 Ui 7 and yi, for i = 1, . . . , n, respectively. Using this notation, the totality of all systems (1) 
is given by 

x = f(x,u,w) - 2 - 
y = h(x,w) 

and set X = W 1 x . . . x R r ". 

The control objective is to reach output agreement of all nodes in the network, independent of the exact 
representation of the time-varying external signals. Therefore, a controller of the form 

Afe = H k (£k), k = 1,2, ...,m, 
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with state G R" fc and input «|; e R p , is placed between any pair of neighboring nodes. When stacked 
together, the controllers (3) give raise to the overall controller 

A = H(Q, W 

where £ G S = R" 1 x . . . x M^™ . Throughout the paper the following interconnection structure between 
the plants, placed on the nodes of Q, and the controllers, placed on the edges of Q is considered. A 
controller (3), associated with edge k connecting nodes has access to the measurement of the relative 
outputs yi — yj . In vector notation, the relative outputs of the systems are 

z = (B®I p ) T y. (5) 

The controllers are then driven by the systems via the interconnection condition 

v = -z, (6) 

where v are the stacked inputs of the controllers. Additionally, the output of a controller A^ influences 
its two incident systems. Thus, the stacked vector of controller's output A drives the process (2) via the 
interconnection 1 

u = (B®I p )\. (7) 
We are now ready to formally introduce the output agreement problem. 



Definition 1 The output agreement problem is solvable for the process (2) under the interconnection 
relations (5), (6), (7), if there exists controllers (4), such that every solution originating from W x X x S 
is bounded and satisfies lim^oo (B T (g> I p ) y(t) = 0. 

The first step is to investigate necessary conditions for the output agreement problem to be solvable. 
The closed-loop system (2), (5), (6), (7), (4) can be written as 

w = s(w) 

x = f(x,(B®I p )H(£),w) (8) 
£ = F(£,~(B®I p ) T h(x,w)). 

If the output agreement problem is solvable, then all the solutions to (8) starting from W x X x S 
converge to the w-limit set Q(W x X x 5). Notice that by boundedncss of the solutions such a set is 
non-empty, compact and invariant. Furthermore, the w-limit set fl(W x X x S) must be a subset of the 
set of all pairs (w, x) for which (B <g> I p ) T h(x,w) — 0. Let now (w,x w ,£ w ) be a solution to (8) starting 
from f2(W x X x £). By the invariance of the w-limit set, we have 



x w = f(x w ,u w ,w) 
= (B®I p ) T h(x w ,w) 



Tui„ w ^ ( 9 ) 



and 

£w = F(f,0) 



u w = (B®I P )H(£ W ). ( 10 ) 



We summarize the necessary condition as follows. 



Proposition 1 If the output agreement problem is solvable, then, for every initial condition in fi(W x 
X x E), there must exist solutions (w,x w ,^ w ) such that the equations (9), (10) are satisfied. 

The constraint (9) ensures that there exists a feed-forward control input u w that keeps the systems in 
output agreement. The second constraint (10) ensures that the controller (4) is able to generate this 
feed-forward input signal. 

1 The interconnection structure (5), (7) naturally represents a canonical structure for distributed control laws. This 
structure is often considered in the context of passivity-based cooperative control, see e.g., [Arc07], [BAW11], [vdSM12], 
[DJ12a], [BZA13a]. 
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The constraints (10) can be rewritten independently of the controller ([IM10]). Let in the following X w 
be some trajectory satisfying u w — (B <E> I q )X w , where u w is a solution to (9). Note that X w is only then 
uniquely defined if the graph Q has no cycles. Otherwise, it can be varied in the circulation space of Q. 
Bearing in mind the structure of the controllers (4), it descends from the constraints (10) that 



e 

A" 



= nr,o) 



(ii) 



Suppose now, that there exists an integer d and maps r : W i-> 
satisfying 

dr 

-s(w) — (j)(r{w)) 



and tp : 



dw 



A" 



ip(r(w)). 



(12) 



Observe that T,</>,ip do not depend on the controller since X w depends on B and u w , the latter being 
dependent only on the process to control and B on the topology of the underlying graph. Now, the 
dynamical system 

fj = <j>(rj), f] £ R d 

a = i>{r,). 



(13) 



has the following property. If r] = t(w(0)), then the solution rj(t) to (13) starting from r] is such that 
X(t) = X w (t) for all t > 0. 

For designing a controller with the structure (4), i.e., that decomposes into controllers on the edges of Q, 
we introduce a vector % £ R d for each edge k = 1, . . . , m, and denote with tpk the entries of the vector 
valued function ip corresponding to the edge k. Each edge is now assigned a controller of the form 



Xk 



ipkiVk), k 



1,2,..., m. 



(14) 



With the stacked vector r] = [r){ , . . . , rj^y and the vector valued functions 





















X^m). 







(15) 



the overall controller is given as 



n = 4>{n) 
x = i,( v ). 



(16) 



Note that if the initial condition is chosen as r]o = I m ® r ( w (0)) then the solution r](t) to (13) starting 
from rjo is such that X(t) = X w (t) for all t>0, which is the same property expressed by (11). 



Remark 1 If the functions 4>(j]) and ip(r]) are linear, that is <f>(r}) — <&rj ($ £ R dxd ), and ip(rj) = ^rj £ 



xd 



), then ip k 



) are the rows of the matrix ^ corresponding to the edge k, and 



the global functions (15) are given by <j>(r)) 
[De 13]). 



(I m ® $)r], and ip(r)) = block.diag(* 



T 
1 ) 



*4) ? 7 ( see a l so 



In summary the overall controller (16) can be interpreted as internal-model-based controllers placed at 
the edges of the graph Q, where all controllers have the same global internal model. The role of the 
internal model in problems of coordination in networked systems has been investigated in [WSA11] for 
linear systems and in [WiclO], Chapter 5 for nonlinear systems. The result above adds up to these 
results. 



Remark 2 The condition (12) might be difficult to meet for general nonlinear problems. However, the 
condition can in fact always be met if both the internal model and the desired feed-forward input are 
linear in the disturbance, i.e., s(w) = Sw and X w = Aw, see [Fra76]. 
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Remark 3 Suppose that the w-limit set can be expressed as 0(W x X xE) = {(w, x,£) : x = n(w), £ — 
tt c (w)}. Then x w = n(w) and the so-called regulator equations (9) express the existence of an invariant 
manifold where the "regulation error" (B T ® I p )y is identically zero provided that the control input u w 
is applied. The conditions (11) express the existence of a controller able to provide u w . Moreover, (9), 
(11) take the familiar expressions, see e.g. [IB90]: 

— s(w) = f(w(w),(B ®I p )H(ir c (w)),w) ^ 
= (B <g> I p ) T h{ir(w),w) 



and 



-^L s {w) = F(n c (w),0) (lg) 
AH = H(tt c (w)). 

Example 1 Consider linear systems of the form 

Xi — AiXi -\- Gilli -\- PiWi ( -\ c\\ 

ii- (It- \ ' 

Define the stacked system with n] T , and A = block. diag{Ai, . . . , An}. Let u, w, G, and 
S, be defined equivalently. There exist x w = Hw and u w — Tw satisfying (9) if and only if II, T satisfy 
the Sylvester equation 



US = AH + GT + P 

(b ® i p fcn = o. 

The equations are feasible if and only if 



TS-IT-T „ (20) 



A - XI G 
is full-row rank for each A in the spectrum of S. 2 



{B ® I p ) T C 



3 Output agreement under time-varying disturbances 

In this section we highlight sufficient conditions that lead to a solution of the problem for a special class 
of systems (1), namely incrementally passive systems, see e.g., [PM08], to which we refer the reader for 
the definition of a regular storage function. 

Definition 2 The system (1) is said to be incrementally passive if there exists a C 1 regular storage 
function V : R>o x R Ti x R n n> K> such that for any two inputs Ui,u\ and any two solutions Xi,x\, 
corresponding to these inputs, the respective outputs yi, y[ satisfy 

dVi dVi dVi , , 

^ + ^fi(x i ,u i ,Wi) + —fi(x i ,u i , W i) (21) 
< ( yi - yi) T (ui - «{)• 

Example 1 (ctnd.) Linear systems of the form (19) that are passive from the input Ui to the output 
satisfy the assumption above, with Vi — \{xi — x' i ) T Qi(x i — x'A and Qi = Qf > the matrix such that 
Af Qi + QiAi < and Qid = Cj . 



2 We can establish at this point a connection to the equilibrium independent passivity framework studied in [BZA13a]. 
In the case A; is invertible and Si = 0, the matrix 11^ that satisfies (20) for any I\ is El; = — A7 1 (BiYi + Pj) and 
xf = -AT 1 (B l T i + Pi)wi. The output y w becomes -C i (AT 1 {B l T i + P x )w x = -C l A~ 1 B i u™ - C l A~ 1 P i w l which 
coincides with the equilibrium input-to-output maps k Vi (ui) considered in [BZA13a], Example 3.4. 
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Example 2 Nonlinear systems of the form 



ii = fi(xi) + dui + PiWi 

Vi — CiXi 

with fi(xi) — VFi(xi), Fi(xi) twice continuously differentiable and concave, and Gi = Cj are incremen- 
tally passive. In fact, by concavity of Fi (a;,), (x i -x-) T (/ i (x i ) — < 0, and Vi = \{xi—x' i ) T (x i -x' i ) 
is the incremental storage function. 

In the previous section, it was shown that the controllers at the edge have to take the form 

m = 4>(m) (99] 

A fe - Mm), fc = l,2,...,m. [ ' 

Now, they must be completed by considering additional control inputs that guarantee the achievement 
of the steady state. While we require the internal model to be identical for all edges, i.e., <j>{r)k), the new 
augmented systems might be different. Then, the controllers modify as 

Vk = <j>k{Vk,Vk) , 23 x 
Afe = ipkiVk), fc = l,2,...,m, 

where all controllers reduce to the common internal model if no external forcing is applied, i.e., <fik{Vk, 0) = 
4>(r)k)- The following is the main standing assumption that the controllers must satisfy to solve the output 
agreement problem for the class of incrementally passive systems: 

Assumption 1 For each k — 1,2, ... ,m, there exists regular functions Wk(r]k,r]' k ), with Wk : x 
R qk ->■ R+ such that 

-g^Mvk,vk) + M Mvk,v k ) (24) 

< (Afe - \' k ) T {v k - v' k ). 

It is in general difficult to design the incrementally passive controllers above. A first simple example 
when the design is possible is when the feedforward control input is linear, that is (12) is satisfied with 
r = Id, 4> = s and ip is a linear function of its argument. In this case, we let 

(t>k(Vk,0) = tpk(Vk) = M k n k 

and define 

<j>k(Vk,Vk) = s(r]k) + M k Vk- (25) 

Then by definition of s as the gradient of a concave function, the storage function Wk(r]k, Vk) — ^("tk — 
VkViVk - v'k) satisfies 

-g^MVk,v k ) + -^rMVk,v k ) = 

(Vk - Vkf (A w k) - s(w k )') (26) 
+(% - v' k ) T Ml(v k - v' k ) 
< (ipk{Vk) - ipk(v'k))( v k - v'k), 

that is (24). 

We state below the main result of the section that, while extending to networked systems the re- 
sults of [PM08], provides a solution to the output agreement problem in the presence of time- varying 
disturbances. It is a slightly more general statement than Theorem 2 in [De 13]. 
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Theorem 1 Consider the system (2) 

w — s(w) 



s (27) 

X = J(X,U,W) 

and let the regulator equations (9) hold. Consider the controllers 

V = Hv,v) , 9R * 

A = ij)(n) + v { ' 

where v is an extra feedback to design, and let (p and ip be the stacked functions of <frk{Vk,Vk) o-nd ipk(Vk) 
with the internal model property being satisfied. Consider the interconnection conditions 



v = -{B T ®I p )y. (29) 



u = (B ®_7 P )A 
If Assumption 1 holds and 

v = —z, 

with z = (B T (8> I p )y, then the output agreement problem is solvable, that is every solution starting from 
W x X x S is bounded and 

lim z(t) = lim (B ® I p ) T y{t) = 0. 

t— >+oo t— >+oo 

Proof: By the incremental passivity property of the x subsystem in (27) and (9), it is true that 

dV dV dV 

<(y-y™f( U -U w ), 

where V = J2i ^i- Similarly by Assumption 1, the system (28) satisfies 

dW - dW - * 

-^-M. v ) + d^^ w ) < (A - \ w ) T v -{v- i/>, 

with W = J2 k Wk and 4>(n w ) — I m <g> 4>(r] w ). Bearing in mind the interconnection constraints u — 
{B®I P )\, u w = (B®I P )\ W , andu = -{B T ®I p )y, and letting U((x, x w ), (rj, n w )) = V(x, x w ) + W(n,n w ) 
we obtain 

U((x, x w ), (V, V w )) ■■= V(x, x w ) + W( V , n w ) 
= (y - y w ) T (u - u w ) + (A - \ w ) T v - (v - v w ) T v 
= {y-y w ) T {B®I P ){\-\ w ) 
- (A - \ W ) T {B T ® I p )y + u w ) T {B T ® I p )y. 

By definition of output agreeement, (B £g) I p ) T y w = and the previous equality becomes 

U{{x,x w ),{fi,rf°)) - v T (B T ®I p )y 

= -\\(B T ®I p )y\\2 = -z T z, 

by definition of v = — z and v w — 0. Since U is non-negative and non-increasing, then U(t) is bounded. 
As x w ,n w are bounded 3 and U is regular, then x, 77 are bounded as well. Hence the solutions exist for 
all t. Integrating the latter inequality we obtain 

r+oo 

/ z T (s)z(s)ds < [7(0). 
Jo 

By Barbalat's lemma, if one proves that j^z T (t)z(t) is bounded then one can conclude that z T (t)z(t) — > 0. 
Now, z(t) = (B T ® I p )y = (B T ® I p )h(x,w) is bounded because x,w are bounded. If h is continuously 



s By definition, (w, x w , rj w ) belongs to the cj-limit set, which is compact. Hence, x w ,r) w are bounded. 
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diffcrcntiable and x,w are bounded, then z is bounded and one can infer that -^z T (t)z(t) is bounded. 
By assumption, w is the solution of w = s(w) starting from a forward invariant compact set. Hence, 
both w and w are bounded. On the other hand, x satisfies 

x = f(x, (B ® I p )i>{ri) - z, w) 

which proves that it is bounded because x, r\, z were proven to be bounded, while w is bounded by as- 
sumption. Therefore, x, w are bounded and this implies that j^z T (t)z{t) is bounded. Then by Barbalat's 
Lemma we have lim t ^ +00 z(t) — as claimed. 



3.1 Linear systems and distribution networks 

We investigate next the output agreement problem for linear dynamical systems and focus on a routing 
control problem in inventory systems under time-varying demand and supply. Consider an inventory 
system with n inventories and m transportation lines, and let B be the incidence matrix of the trans- 
portation network. The dynamics of the inventory system is given as 

x = BX + Pw, (30) 

where x € K n represents the storage level, A 6 M m the flow along one line, and Pw an external in- 
/outflow of the inventories, i.e., the supply or demand. We assume that the time varying supply/demand 
is generated by a linear dynamics 

w = Sw (31) 

and that it is balanced at any time, i.e., l T Pw(t) — for all t > 0. The control problem is to find a 
distributed control law of the form 

m = ^kVk + AfeUfc 
Afe = *fc% + TfcUfc 

such that for each initial condition (wo, xq, 770) the solution of the closed loop system remains bounded 
and a balancing of the inventory levels is achieved, i.e., lim t _ i . 00 B T x = 0. 

The closed-loop system (30), (31), (32) can be understood as feedback interconnection (5), (6), (7), 
of the controller (32) and the linear system 

w = Sw , % 

• P (33) 

X = u + Pw, 

i.e., with A = and G = I n . Thus, the distribution problem can be understood as an output agreement 
problem. The regulator equations (9) become 

US = T + P, = B T Il. (34) 

These conditions are satisfied with II = J, i.e. the all-ones-matrix of appropriate dimension, and T = 
JS — P. We can directly conclude from Remark 2, that there exists a matrix H such that X w = Hw. 4 
The steady-state solution is such that x w = Hw = Jw = x™ 1. Since, by assumption, l T Pw = at any 
time instant, it holds that l T x(t) = for all t>0. However, this implies that l T x w = $1 T 1 — 0. The 
latter condition can only be satisfied for $ = 0, and therefore the steady-state solution x w must be a 
constant vector. Now, the steady-state routing X w satisfies 

= BX W + Pw. (35) 

Thus, any controller (32) solving the output agreement problem, solves at the same time the exact 
routing problem of the time- varying supply /demand. The condition (12) is satisfied with t = Id, <fi = S 
and \p = H. Following our previous discussion, it remains to design the incrementally passive local 



4 How such a matrix H can be found is discussed in ([De 13]), and we refer the interested reader to this reference. 
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controllers. Let in the following H £ denote the fc-th row of the matrix H . We know from (25), that the 
internal model controller on the edge k of the form 



Vk = Srj k + H k v k 
= H k n k 



(36) 



is incrementally passive. Finally, it follows directly from Theorem 1 that the distributed internal model 
controller 



f) = Sr)- HB T x 
X = Hr]- B t x 



(37) 



where S — I n <g> S and H = block. diag(i/j r , . . . ,-ff„), solves the output agreement problem and, addi- 
tionally, achieves an exact routing of the time- varying supply through the network. This result has been 
also proven in [De 13], but is presented here in the more general context of output agreement. 

3.2 Output agreement in the case of constant disturbances 

The critical assumption in the derivation presented above seems the incremental passivity of the internal 
model, i.e., Assumption 1. However, the proof above exploits Assumption 1 only in a weaker form. 
In particular, it requires the incremental passivity property (24) not to hold with respect to any two 
trajectories, but only with respect to the real and the steady-state trajectory, i.e., with n' k = rfg,v' k — 
0, X' k — X w . Bearing in mind this observation, it is possible to find a storage function W k that fulfills (24) 
in the case of constant disturbances. In fact in this case (r] k ,X k ) satisfy X k = tpk(Vk) f° r some constant 
r, w , i.e., rfg = 0. 

Let now ip k be a strongly monotone function, and consider the following storage function ([JOGC07], 
[BZA13a]): 



W fe ( % ,0=*fe(%)-* fe (0 

+w£W)(%-^), 



where : R 9 — >■ R is a twice continuously differentiable function such that V^fkiVk) = tpk(Vk)- 5 Now 
if ipk is monotone, <3/fc is convex and, by the global under-estimator property of the gradient [BV03], we 
have 

for each r\ kl rf^. If ^ k is strictly convex, i.e., ip k is strongly monotone, then the previous inequality holds 
if and only if r/k — f] k ([BV03]), and W k is regular ([JOGC07]). Hence, W k is a non-negative storage 
function if ^ k is convex, and a positive regular storage function if is strictly convex. Furthermore, 



dr) k 



(f>k(Vk,v k ) = (ipk(Vk) ~ V'fe(O) v k- 



Hence, in the case of constant disturbances, Assumption 1 is always fulfilled and the output agreement 
problem is solved by the controllers 



f) = v_=-(B®I p )y 
X = $( V ) - (B T ® I p )y. 



(39) 



Control laws of the form (39) have been studied in the context of network clustering in [BZA13b] , [BZA11] 
and in a more general network optimization framework in [BZA13a]. 

4 Time-varying optimal distribution problems 

We revisit now the distribution problem of inventory systems discussed in Section 3.1, i.e., 

x = BX + Pw, (40) 



3 Note that W k is the Bregman distance between r) and r) w for the function [Brc67]. 
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with a time- varying external demand/supply. Let for now the supply/demand vectors be generated by 
a possibly nonlinear dynamics w = s(w). In contrast to Section 3.1, the control objective is not only to 
balance the inventory levels, but additionally to achieve an optimal routing in the network. We associate 
to each transportation line a cost 

Pfc(Afc), fc = l,2,...,m 

for the flow Afe. Let Vk be convex and continuously differentiable. We aim here to design routing 
controllers on the transportation lines of the form (22), taking only the imbalance between the incident 
inventories as inputs such that asymptotically all inventory levels are balanced and the flows in the 
network minimize the flow cost defined by the cost functions Vk- 

Before moving to the dynamic control problem, we briefly review the static optimal distribution problem. 
Consider a fixed constant supply vector w. The (static) optimal distribution problem is to find a routing 
X w G R m such that 

m 

X w = argmin ^"Pfc(A fe ) 

k=i (41) 

= BX + Pw. 

In the following, the notation V(X) = Y^k=i ^M^fc) wm be used. The Lagrangian function associated to 
(41) with the multiplier £ is 

C(X, v) = V(X) + ( T (-BX - Pw). 

From the Lagrangian, one obtains directly the optimality conditions (KKT conditions). In particular, 
(X w ,v w ) is an optimal primal/dual solution pair to (41) if the following nonlinear equations hold 

VV(X W ) - B T C =0 

v ; (42) 
BX W +Pw = 0. J 

Note that the first condition simply implies W(X W ) G 1Z(B T ) since £ w has no further constraints. Wc 
define the optimal routing/supply pairs as 

r = {(X,w) el m xW : VP(A) G K(B T ), 

BX + Pw = 0}. 

In particular, (X w ,w) G T if and only if X w is an optimal solution to the static optimal distribution 
problem (41) with the supply vector w. The main difficulty associated with the set T relates to the 
constraint VP(A) £ 1Z(B T ). This constraint can be avoided if the optimality conditions are expressed 
in terms of the dual solutions £. In what follows, we impose the condition that VP is invertible. 6 The 
two optimality conditions (42) can now be expressed as the following single nonlinear expression 

BVV- 1 (B T ( w ) + Pw = 0. 

Bearing this in mind, we define the set of optimal dual solutions as = {(£, w) G W l x W : 
BVV^ 1 (B T () + Pw = 0}. We want to emphasize two properties of To- First, if (( w ,w) G then 
(C w +cl,w) G Yd for any c G M. Second, (C t w , w) G if any only if the corresponding routing strategy 
X w = VV^ 1 (B T ( W ) satisfies (X w , w) G T. We are now ready to formalize the dynamic problem. 

Definition 3 The time- varying optimal distribution problem is solvable for the system (30), if there 
exists a controller (22) such that any solution originating from W x X x S is bounded and satisfies (i) 
linit^co B T x(t) = 0; and (ii) lim^oo distr(A(i), w(t)) = 0. 

Our previous discussion revealed that it can be advantageous to express the optimality conditions in 
terms of the dual solutions £. In fact, we use this observation here to refine the control and restrict our 
attention to dynamic controllers of the form 



X = \>p- 1 (B T iP(r 1 )). 



(43) 



3 Note that this is always given if V is strongly convex. 
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This control structure can be interpreted as an internal model controller for the dual variables of the 
optimal network flow problem. For the time-varying optimal distribution problem to be feasible, it is 
necessary that the manifold 



H(w,x,r]) = 



B T x 

BW- 1 {B T i;(T]))+Pw 



= (44) 



is invariant under the closed-loop dynamics 

w = s(w) 

x = BVV- 1 (B T i)(r]))) + Pw (45) 
rj = <p(rj, B T x). 

Note that, in contrast to the original output regulation problem, the "output" function H depends 
explicitly on the state of the controller r). However, at this point, the advantage of the internal model 
controller design for the dual variables becomes obvious. Let (w, x w , r} w ) be a solution to (45) starting 
in f2(W x X x S), satisfying 

h(x) -~ B T x w = 0. (46) 

The previous condition implies x w = j31. Now, by the structure of the inventory dynamics follows 
l T x = at any time, and consequently 

x w = = BVP- 1 (B T ^( V W ))) + Pw. 

Thus, the corresponding routing strategy X w — \7P~ 1 (B T ip(rj w ))) is optimal at any point in time, 
i.e., (X w (t),w(t)) G T. Thus, by restricting the discussion to the "dual" controller structure (43), we 
transformed the time- varying optimal distribution problem into an output agreement problem. 
If the time- varying optimal distribution problem is solvable, then for any solution of the system w = s(w) 
there exists a solution to ij w satisfying fj w = (j){ , q w ) such that 

BVV- 1 {B T ^(r ] w '(*)))) + Pw(t) = 0. (47) 

Under sufficient smoothness conditions on P^ 1 and on s(x), there exists a differentiable trajectory rf 
satisfying (47). Keeping this observation in mind, we can formalize conditions on the controller (43). 
Suppose there exists a map t:W4H such that 

dr 

— s{w) - 0(r(«;)), Vw e W, 
BVV- 1 (B t iP(t(w))) +Pw = 0. 

Then, the dynamical system i) = (/)(r)), \i = ip(ri) has the following properties. If 770 = t(w(0)), then the 
solution rj(t) starting at rj is such that tp{rj{t)) = ( w 7 that is (^(t),w(t)) e T D . Thus, if the condition 
(48) is met, the internal model controller generates the correct time-varying dual solutions, and can 
produce the optimal routing strategy. We formalize this in the following result. 



Proposition 2 Let condition (48) hold. Then, for any initial condition w(0) € W, there exist an initial 
condition (xo,r]o) G X x S such that the solutions (w(t),x(t),n(t)) to (45) satisfy (i) B T x(t) — and 
(11) (X(t),w(t)) G T 7 where X(t) = \7p- 1 (B T i;(r 1 (t))). 

Condition (48) is the counterpart to the original controller regulator equation (12). Note that, while 
the regulator equations for the systems dynamics are easily satisfied, as we discussed in Section 3.1, the 
difficulty of the time- varying optimal distribution problem results from partially fixed structure of the 
controller (43). We discuss next the distributed controller design for two specific representations of the 
optimal distribution problem. 
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4.1 The Linear-Quadratic Case 

Suppose the supply is generated by a linear system 

w = Sw, S + S T = 0, 
and the transportation cost functions are quadratic 

= ^ T QA, 

for Q = diag(<7i, . . . , q m ) and g fe > 0. In the linear-quadratic case one can solve the conditions (48) 
directly. Consider the conditions (48) with t(w) = w, and iP(t{w)) = Hw, i.e., 

BQ- 1 BHw + Pw = 

q , (49) 

Sw = <p{w). 

This requires <p(rj) = Sr/, i.e., the controller is a copy of the exosystem. It remains to find a matrix H. 
Note that Lq — BQ~ 1 B T is a weighted Laplacian matrix of the network. As Lq has one eigenvalue at 
zero, with the corresponding eigenvector 1, it is not invertible. However, one possible solution to (49) is 

H = -L Q P, 

where Lq pseudo-inverse of the weighted Laplacian Lq, see e.g., [GX04]. 7 We can now construct 
a distributed controller. Recall that (43) provides an internal model controller for the dual variables, 
assigned to nodes of the network. Thus, we assign an internal model controller to each node and introduce 
the variable rji € R d , satisfying the dynamics 



(i = Hjm, i = l,...,n 



where Hf is the i-th row of H. The routing at one edge computes then simply as Afc = q^iCj — Ci ) 7 
where j, i are the nodes incident to edge k. After introducing 77 and £ as the stacked vectors of rji and 
Ci, respectively, and the matrices S — (/„ <g> S), H = block. diag(_ffj r , . . . , H^), we can express the overall 
controller as 

fj = Sr) 

X = Q~ 1 B T ( = Q~ 1 B T Hn. 

This distributed internal model controller needs to be augmented with additional control inputs to ensure 
convergence to the desired stead state. This is formalized in the next result. 

Theorem 2 Consider the inventory system (30) with the supply generated by the linear dynamics w — 
Sw. Consider the controller 

77 = Srj + H T BQ~ x v 
X = Q- 1 B T Hr 1 + v. 

with the interconnection condition v = —B T x. Then, every solution of the closed-loop system is bounded 
and (i) lim t ^ +oc B T x = 0, and (ii) lim 4 _j. +00 distr(A(f), w(t)) = 0. 

The proof of this result follows directly along the same lines as the proof of Theorem 1, taking into 
account that V(x 7 x w ) = \\\x — x w \\ 2 and W — \\\r] — n w \\ 2 are regular incremental storage functions. 
For concluding the optimality, one has to note that any trajectory in the set £ = {{x,rf) : B T x = 0}, 
which is invariant under the dynamics, is such that x — 0. Thus, any trajectory fj in £ is such that the 
corresponding routing 

A = Q~ 1 B T Hr) = VV- 1 {B t Ht)) 

satisfies BX + Pw = 0, i.e., the optimality conditions (42). Thus, any routing A in the invariant set £ is 
an optimal routing. 

7 By the properties of L^q, it is promptly verified that 

BHw + Pw = -BQ- 1 B T Lt ) Pw + Pw = 

-(/ )Pw + Pw = 0. 

n 
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4.2 Optimal distribution with constant supply 

A second version of the optimal distribution problem, that can be solved by the internal model approach, 
relates to problems with constant supply and demand, i.e., s(w) = 0. In this case, we need no restriction 
to quadratic objective functions, but can consider general convex cost functions Vk- The conditions (48) 
become 

= 0(t(w)) 
BVV- 1 (B t iP(t(w))) +Pw = 0, 

and they are solved with <j>(w) — and t(w) such that (t(w),w) G Tp. In the static case, it becomes 
again advantageous to place the internal model controllers at the edges, and not at the network nodes. 
However, we still consider an internal model for the dual variables, but consider now the dual variables 
a = B T ( instead of £ as introduced in (42). Consider now the controller 

a = u, cr(0) G Tl{B T ) 
\ = \7T- 1 (a) + v. 

Note that if v{t) G K(B T ) then also a(t) G TZ(B T ) for all times t > 0. If the initial condition a(0) is 
chosen as a w — B T ( W , with (( w ,w) G Td, then, under the input v = 0, the controller (52) generates the 
desired output. The internal model controller (52) is now not incrementally passive for arbitrary input 
signals, but it is incrementally passive with respect to any constant input signal. A storage function can 
therefore be found in the structure (38). 8 Based on our previous discussions, we can now immediately 
conclude that the controller (52) with the interconnection condition v = —B T x solves the optimal 
distribution problem. 

5 Conclusions 

We proposed an internal model control design approach for output agreement of incrementally passive 
nonlinear systems with time- varying external disturbances. Building upon theses results, we studied the 
optimal distribution problem in inventory systems with time- varying supply and demand. These prob- 
lems can be cast as output agreement problems if their dual formulation is considered. We showed how 
the optimal distribution problem can be solved by internal model controllers for the dual variables. The 
specific solution to the distribution problem is discussed for the linear-quadratic case and for problems 
with constant supplies. 

In the case of constant disturbances and equilibrium independent passive systems, the paper [BZA13a] 
characterizes the value on which the outputs agree and establishes insightful relations with a number 
of network optimization problems. In the case of time-varying disturbances, such a characterization is 
still elusive and represents an interesting future research avenue. In the same work, the role of bounded 
interactions is investigated, a topic that has not been addressed in this paper. Passivity has provided a 
natural setting to study cooperative control algorithms in the presence of delays. The effect of delays in 
problems such as those studied in this paper has been so far neglected, to the best of our knowledge. It 
would be interesting to fill in this gap. 
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